#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <map>
#include <algorithm>
#include <fstream>
#include <time.h>

#include "fastjet/ClusterSequence.hh"
#include "fastjet/tools/Pruner.hh"

#include "Njettiness.hh"

#include "CmdLine.hh"
#include "TFile.h"
#include "TChain.h"

#include "TH1F.h"
#include "TH2F.h"

#include "Analysis.h"


class FileSkimmer : public Analysis
{
public:

  TTree* treeout;
  FileSkimmer(const CmdLine& cmdline)
    : Analysis(cmdline)
  {
    fileout->cd();
    treeout=new TTree("STDHEP", "STDHEP");
    hist["met"]=new TH1D("MEt", "#slash{E}_{T}", 50, 0, 500);
    //    hist["nevt"]=new TH1D("nevt", "nevt", 10, 0, 10);
    //    hist["nevt_unweighted"]=new TH1D("nevt_unweighted", "nevt_unweighted",
    //				     10, 0, 10);


    treeout->Branch("weight", &weight, "weight/D");  
    treeout->Branch("size", &size, "size/I");
    treeout->Branch("PID", PID, "PID[size]/I");
    treeout->Branch("Status", Status, "Status[size]/I");
    treeout->Branch("E", E, "E[size]/D");
    treeout->Branch("Px", Px, "Px[size]/D");
    treeout->Branch("Py", Py, "Py[size]/D");
    treeout->Branch("Pz", Pz, "Pz[size]/D");
    treeout->Branch("Pt", Pt, "Pt[size]/D");
    treeout->Branch("Eta", Eta, "Eta[size]/D");
    treeout->Branch("Phi", Phi, "Phi[size]/D");
  }

  bool cut()
  {
    hist["nevt"]->Fill(0.0,weight);
    hist["nevt_unweighted"]->Fill(0.0);


    if(cmdline.present("-Mt"))
      {
	float Mt_cut=cmdline.value<float>("-Mt", 400);
	float Mt_present=0;
	for(int i=0; i<jets.size(); i++)
	  {Mt_present+=jets[i].m();}
	if(Mt_present < Mt_cut)
	  return false;
      }

    if(cmdline.present("-Pt"))
      {
	float Pt_cut=cmdline.value<float>("-Pt", 200);
	if(jets.size()<1)
	  return false;
	if(jets[0].pt() < Pt_cut)
	  return false;
      }


    if(cmdline.present("-jet4"))
      {
	if(Event_==0)
	  {
	    cout<<"OPTIONS: Demand 1j > 100, at least 4j > 50"<<endl;
	    hist["nevt"]->GetXaxis()->SetBinLabel(2, "jet4");
	  }
	
	if(jets.size()<4)
	  return false;
	if(jets[0].pt()<100)
	  return false;
	if(jets[3].pt() <50.0)
	  return false;
      }

    
    hist["nevt"]->Fill(1.0,weight);
    hist["nevt_unweighted"]->Fill(1.0);
    return true;
  }

  void analyze();

};



int main(int argc, char** argv)
{

  //initialize command line interface
  CmdLine cmdline(argc,argv);

  if(!cmdline.present("-i") || !cmdline.present("-o"))
    {
      cout<<"FileSkimmer: trim down rootfile by applying predefined cuts";
      cout<<"Usage: -i <input> -o <output> -a [algorithm] -r [radius] -m [mass] -n [nev] -p -w"<<endl;
      cout<<"-a \t clustering algorithm to be used"<<endl;
      cout<<"-r \t clustering radius"<<endl;
      cout<<"-n \t number of events to process"<<endl;
      cout<<"-p \t whether jet pruning is used"<<endl;
      cout<<"-w \t whether event has a weight stored as evt Number"<<endl;
      cout<<"-gluino \t retrieve plots that trace the gluinos "<<endl;
      cout<<"-jet4 \t apply 4 jets cut "<<endl;
      cout<<"-jet2 \t apply 2 jets cut "<<endl;
      cout<<"-display \t number of display plots"<<endl;
      
      return 1;
    }

  //create analysis object
  FileSkimmer myanalysis(cmdline);
  myanalysis.tree_init();
  myanalysis.analyze();

  return 0;
}


//main function to loop over all the entries
void FileSkimmer::analyze()
{


  //begin event loop
  for(int entry=0; 
      (maxevt<0 || entry<maxevt) && entry<tree.GetEntries();
      entry++)
    {

      if(entry%1000==0 && entry/1000 > 0)
	cout<<"MESSAGE: "<<entry<<" Processed"<<endl;

      //populate our variables
      tree.GetEntry(entry);

      Event_=entry;
      if(cmdline.present("-unit-weight"))
	{
	  weight=1.0;
	}

      makejets();
      
      //grab all the jets
      if(!cut()) continue;
      treeout->Fill();
      
    }
  
  treeout->Write();
  fileout->Write();
  fileout->Close();
}


